Hedonic
  • Home
  • Risk
  • Amenities
  • Trails
  • Beaches

On this page

  • OSM Beaches
    • Matching Polygons
  • Coral Reefs
    • Matching Polygons

  • Show All Code
  • Hide All Code

  • View Source

OSM Beaches

We are going to use a beach shapefile from OpenStreetMaps and calculate the length and the width of the identified beaches within the maps. We will also need to link these to Hawaii poylgons like what was done in the buffer approach. Here what we do is use a buffer to link houses that adjoin the polygon beaches

Code
beaches <- read_sf("/Users/ashleylowe/Trails/beaches_meters.shp")
beaches <- st_transform(beaches, 32605)
beaches_buffered <- st_buffer(beaches, dist = 10)
beaches_buffered <- st_transform(beaches_buffered, 4326)
source("_common.R")



intersections <- st_intersects(HI_2023_shapefile, beaches_buffered)
HI_intersections <- HI_2023_shapefile[lengths(intersections) > 0, ]

tmap_mode("view")

center_lon <- -157.45
center_lat <- 20.75
zoom_level <- 7

tm_shape(beaches) +
    tm_polygons(
    col = "name",
    palette = "orange",
    size = 0.1,
    id = "Name",
    popup.vars = c("Name" = "name"),
    legend.show = FALSE      # hides the legend
  ) +
  tm_basemap() +
  tm_view(set.view = c(center_lon, center_lat, zoom_level))

Matching Polygons

Using the Adjoining properties of these beaches we will now merge the beach names in OSRM (note not all beaches have a name), beach width and length based on the polygon calculate using the skeleton of the polygon in QGIS.

Code
# Create a vector of indices for the first overlapping shp2 feature per shp1 row
HI_2023_shapefile$beachname <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$name[idx], collapse = ", ")
  }
})

# Create a vector of indices for the first overlapping shp2 feature per shp1 row
HI_2023_shapefile$BeachWidth <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$BeachWidth[idx], collapse = ", ")
  }
})

HI_2023_shapefile$Beachlength <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$LENGTH[idx], collapse = ", ")
  }
})

Coral Reefs

I am also interested in the coral reefs and if those are incorporated into housing prices. Using the Asner et al Using the property file we isolate the potential coastal properties following

  • Selecting the coastline shapefile

  • Create a buffer of 10meters for coastline

  • Select Hawaii properties intersecting the coastline

  • Buffer properties 300 meters to grab the areas coral cover average.

The following file shows properties which returned a value using this method.

Code
asner <- read_csv("/Users/ashleylowe/Hawaii - Recreation/Island_Condition/Coral/asner_coastline50meters.csv")

asner$costline_property=1

asner=HI_2023_shapefile %>%
    semi_join(asner, by = "objectid") %>%       # keep only matched rows
    left_join(asner %>% select(objectid, coral_asner), by = "objectid")  # bring in coral_cover
# Transform to WGS84
asner <- st_transform(asner, 4326)

tmap_mode("view")
center_lon <- -157.45
center_lat <- 20.75
zoom_level <- 7

tm_shape(asner) +
  tm_polygons(
    "coral_asner",
    palette = "viridis",
    title = "Coral Reef Cover",
    popup.vars = c("Coral Reef Cover" = "coral_asner")
  ) +
  tm_basemap() +
  tm_view(set.view = c(center_lon, center_lat, zoom_level))

Matching Polygons

Matching these to the housing locations across the islands

Code
# Make sure both have the same CRS
st_crs(PB_HawaiiInfo_joined) == st_crs(asner)  

# If FALSE, transform one:
#PB_HawaiiInfo_joined <- st_transform(PB_HawaiiInfo_joined, st_crs(Special_Management_Areas))

# Spatial join: attach polygon attributes to points
PB_HawaiiInfo_joined<- st_join(PB_HawaiiInfo_joined, asner, join = st_within)
Source Code
---
format:
  html:
    code-fold: true        # Enables dropdown for code
    code-tools: true       # (Optional) Adds buttons like "Show Code"
    code-summary: "Show code"  # (Optional) Custom label for dropdown
    toc: true
    toc-location: left
    page-layout: full
editor: visual
---
```{r include=FALSE}
source("_common.R")
```

# OSM Beaches

We are going to use a beach shapefile from OpenStreetMaps and calculate the length and the width of the identified beaches within the maps. We will also need to link these to Hawaii poylgons like what was done in the buffer approach. Here what we do is use a buffer to link houses that adjoin the polygon beaches

```{r cache=TRUE, echo=TRUE, warning=FALSE, message=FALSE}
beaches <- read_sf("/Users/ashleylowe/Trails/beaches_meters.shp")
beaches <- st_transform(beaches, 32605)
beaches_buffered <- st_buffer(beaches, dist = 10)
beaches_buffered <- st_transform(beaches_buffered, 4326)
source("_common.R")



intersections <- st_intersects(HI_2023_shapefile, beaches_buffered)
HI_intersections <- HI_2023_shapefile[lengths(intersections) > 0, ]

tmap_mode("view")

center_lon <- -157.45
center_lat <- 20.75
zoom_level <- 7

tm_shape(beaches) +
    tm_polygons(
    col = "name",
    palette = "orange",
    size = 0.1,
    id = "Name",
    popup.vars = c("Name" = "name"),
    legend.show = FALSE      # hides the legend
  ) +
  tm_basemap() +
  tm_view(set.view = c(center_lon, center_lat, zoom_level))

```

## Matching Polygons

Using the Adjoining properties of these beaches we will now merge the beach names in OSRM (note not all beaches have a name), beach width and length based on the polygon calculate using the skeleton of the polygon in QGIS.

```{r eval=FALSE}


# Create a vector of indices for the first overlapping shp2 feature per shp1 row
HI_2023_shapefile$beachname <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$name[idx], collapse = ", ")
  }
})

# Create a vector of indices for the first overlapping shp2 feature per shp1 row
HI_2023_shapefile$BeachWidth <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$BeachWidth[idx], collapse = ", ")
  }
})

HI_2023_shapefile$Beachlength <- sapply(intersections, function(idx) {
  if (length(idx) == 0) {
    return(NA_character_)
  } else {
    paste(beaches_buffered$LENGTH[idx], collapse = ", ")
  }
})
```

# Coral Reefs

I am also interested in the coral reefs and if those are incorporated into housing prices. Using the ![Asner et al](<https://www.pnas.org/doi/10.1073/pnas.2017628117>) Using the property file we isolate the potential coastal properties following

-   Selecting the coastline shapefile

-   Create a buffer of 10meters for coastline

-   Select Hawaii properties intersecting the coastline

-   Buffer properties 300 meters to grab the areas coral cover average.

The following file shows properties which returned a value using this method.

```{r cache=TRUE, echo=TRUE, warning=FALSE, message=FALSE}

asner <- read_csv("/Users/ashleylowe/Hawaii - Recreation/Island_Condition/Coral/asner_coastline50meters.csv")

asner$costline_property=1

asner=HI_2023_shapefile %>%
    semi_join(asner, by = "objectid") %>%       # keep only matched rows
    left_join(asner %>% select(objectid, coral_asner), by = "objectid")  # bring in coral_cover
# Transform to WGS84
asner <- st_transform(asner, 4326)

tmap_mode("view")
center_lon <- -157.45
center_lat <- 20.75
zoom_level <- 7

tm_shape(asner) +
  tm_polygons(
    "coral_asner",
    palette = "viridis",
    title = "Coral Reef Cover",
    popup.vars = c("Coral Reef Cover" = "coral_asner")
  ) +
  tm_basemap() +
  tm_view(set.view = c(center_lon, center_lat, zoom_level))


```

## Matching Polygons
Matching these to the housing locations across the islands

```{r eval=FALSE}

# Make sure both have the same CRS
st_crs(PB_HawaiiInfo_joined) == st_crs(asner)  

# If FALSE, transform one:
#PB_HawaiiInfo_joined <- st_transform(PB_HawaiiInfo_joined, st_crs(Special_Management_Areas))

# Spatial join: attach polygon attributes to points
PB_HawaiiInfo_joined<- st_join(PB_HawaiiInfo_joined, asner, join = st_within)

```